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METODO DE VISUALIZACION DE FLUJOS MULTIFASICOS USANDO 
TOMOGRAFIA DE CAPACITANC1A ELECTRICA 



DESCRIPCION 

CAMPO TECNICO DE LA INVENCION 

La presente invention se relaciona con un nuevo metodo de reconstruction 
de imagenes para la visualization de flujos multifasicos usando tomografia de 
capacitancia electrica (TCE), para empiearse en medidores tomograficos de flujo 
multifasico en la cuantificacion de los diferentes fluidos (como gas, aceite y agua) 
producidos por un pozo petrolffero. Tambien se usa para la visualization y 
monitoreo de otros flujos y procesos multifasicos que ocunren en la industria tales 
como: optimization del diseno y operation de separadores, optimization del 
diseno de unidades de desintegracion de fluidizacion catalitica (FCC) y sistemas 
de lecho fluidizado, y optimization del diseno de reactores de lecho fijo. 

ANTECEDENTES DE LA INVENCION 

En terminos generates, la tomografia sirve para obtener una imagen del 
corte transversal de un objeto en un piano determinado. La tomografia de rayos X 
fue la primera en ser desarrollada (en la decada de 1970) y actualmente se usa 
de manera rutinaria en el area medica, asi como en algunas aplicaciones 
industrials (por ejemplo, inspection interna de piezas y detection de fallas en 
materiales). 

Posteriormente, se han venido desarrollando nuevos metodos de tomografia 
orientados mas bien a procesos industrials, conocidos con el nombre generico 
de tomografia de procesos (Williams y Beck, 1995). La finalidad de estos 
metodos, que comenzaron a evolucionar a mediados de la decada de 1980, es 
obtener una imagen de la distribution espacial de las fases o componentes en un 
proceso industrial, usando Onicamente sensores externos y sin causar ninguna 
perturbation en el proceso, tal como se indica en la figura 1 , que muestra un 
sistema de tomografia de procesos compuesto de un tanque o tuberia (A), un 
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sistema de adquisicion de datos (B) y una computadora de reconstruction de 
irnagenes (C). En otras palabras, la tomografia de procesos proporciona una 
manera de , mirar J desde afuera al interior de los mismos, sin necesidad de 
alterarlos fisicamente, y esto representa un modo radicalmente distinto de obtener 
information estructural del proceso a nivel global, a diferentia de los metodos 
tradicionales, que se basan en el muestreo local en un determinado numero de 
puntos. El proceso puede ocurrir en un reactor o mezcladora, en un lecho 
fluidizado, en el interior de un separador, o dentro de una tuberia transportando 
flujo multifasico, por poner algunos ejemplos. 

Existe toda una gama de printipios y tecnicas que pueden seir explotados en 
la tomografia de procesos, incluyendo metodos electricos, basados en la medition 
de impedancia, ultrasonicos, de resonancia magnetica, opticos, y a base de 
radiation ionizante (rayos X o gama). En general, los metodos de radiation 
ionizante producen las im&genes de mejor definition, pero son relativamente 
lentos. Por otro lado, los metodos electricos producen irnagenes de baja 
resolution pero son mucho mas rapidos, de construction robusta y de costo 
relativamente bajo. De ahi la necesidad de contar con mejores metodos de 
reconstruction de irnagenes que, tales como los de la presente invention, 
permitan obtener irnagenes con una mayor definition y fidelidad. 

Particularmente en el caso de la tomografia de impedancia electrica o, 
simplemente, tomografia electrica, se ha dado un progreso muy notable en los 
ultimos afios. Este tipo de tomografia tiene dos modalidades principales: la 
tomografia de capacitancia y la tomografia de resistencia. En un sistema de 
tomografia de capacitancia (Beck et al., 1997; Gamio, 1997; Plaskowski et a/., 
1995), (tal como se indica en la figura 2, formado por un sensor (A), un sistema de 
adquisicion de datos (B) y una computadora de reconstruction de irnagenes (C)), 
normalmente usado con mezclas donde la fase continua es aislante, se utiliza un 
sensor, tal como se indican en las figuras 3 y 4, formado por un arreglo circular de 
electrodes (B) distribuidos alrededor de la section transversal a examinar, y se 
mide la capacitancia electrica de todos los posibles pares de electrodos. Por 
medio de una computadora y un algoritmo de reconstruction de irnagenes 
apropiado, los datos medidos son usados para producir un mapa (o sea, una 
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imagen) de las variaciones de la constante dielectrica (o permitividad relativa) 
dentro del sensor, proporcionando asf una indication grafica de la distribucion 
interna de los diferentes componentes de la mezcla. Con relation a las figuras 3 y 
4, los electrodos (B) pueden estar colocados en el exterior de un tubo aislante (A), 

5 a fin de simplificar la fabrication y evitar el contacto directo con los fluidos (E) del 
proceso. Un segundo tubo metalico externo electricamente aterrizado (C) sirve 
como blindaje o pantalla electrica y para proporcionar resistencia mecanica. El 
sensor tambien cuenta con dos electrodos de guarda cilmdricos (D) que se 
mantienen asimismo aterrizados. 

10 La patente britanica GB 2 214 640, del 6 de septiembre de 1989, describe un 

sistema de tomografia de capacitancia electrica en el cual el metodo empleado 
para realizar la reconstruction de imagenes es un algoritmo de proyeccion inversa 
lineal, conocido en ingles como linear back-projection' (o 'LBP'). Sin embargo, 
dicho metodo de reconstruction produce imagenes de calidad relativamente baja. 

is Los metodos de reconstruction de imagenes propuestos en la presente invention 
permiten superar esta seria limitation. 

La tomografia de resistencia, por otro lado, esta orientada a mezclas donde 
la fase continua es conductora (Plaskowski et a/., 1995; Williams y Beck, 1995). 
En este caso, los electrodos se colocan al ras de la superficie interior del 

20 recipiente que contiene la mezcla, y en contacto directo con los fluidos de la 
misma. Se aplican una serie de corrientes de excitation y se miden los voltajes 
producidos, y a partir de ellos se encuentra la distribucion de conductividad dentro 
del sensor, la cual indica la distribution de los diferentes componentes. 

En principio, la tomografia de capacitancia electrica (TCE) tiene importantes 

25 aplicaciones en la medicion de flujo multifasico, y en particular de flujo bifasico 
gas-aceite, tal como se presenta comunmente a la salida de muchos pozos 
petrolfferos. La manera traditional de cuantificar los diferentes fluidos producidos 
por un pozo consiste en separar la mezcla de los mismos por gravedad usando 
grandes tanques, para posteriormente medir cada componente en forma 

30 individual usando medidores de flujo homogeneo convencionales. En las ultimas 
dos decadas han surgido medidores de flujo multifasico que permiten cuantificar 
la production de un pozo sin necesidad de separar la mezcla (Thorn et ah, 1997). 
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Sin embargo, los medidores multifasicos disponibles actualmente sufren de una 
sensibilidad a cambios en el patron de flujo, a menos que se utilicen con 
dispositivos mezcladores o acondicionadores de flujo, lo cual introduce perdidas 
irrecuperables de presion (que finalmente se traducen en perdidas de energia). La 
limitation anterior se podria superar mediante la TCE, ya que esta permrte 
conocer el patron de flujo y ese conocimiento se podria utilizar para compensar la 
respuesta de los medidores multifasicos convencionales, o, alternativamente, 
tambien es posible en principio disenar un nuevo tipo de medidor multitask*) 
tomografico, basado en el analisis de las imagenes de TCE obtenidas en dos 
pianos transversales de la tuberia proximos entre sf (Hammer ef a/., 1997; 
Plaskowski ef a/., 1995). Adicionalmente, la TCE tiene aplicaciones potenciales en 
la visualization, monitoreo y control de numerosos procesos multifasicos 
industrials. 

Sin embargo, el factor principal que hasta ahora ha limitado la aplicacion 
practica de la TCE es la poca fidelidad o exactitud de las imagenes obtenidas 
usando los metodos de reconstruccion de imagenes disponibles, de ahi la 
importancia de contar con metodos mejorados como el de la presente invention 
(Yang y Peng, 2003). Los metodos directos simples como el de proyection 
inversa lineal (LBP por sus siglas en ingles) producen imagenes relativamente 
deficientes que solo proporcionan una indication cualitativa de la distribution de 
componentes dentro del sensor. Tal es el caso de los metodos que se describen 
en la patente britanica GB 2214640 del 6 de septiembre de 1989. 

Por otro lado, los metodos mas sofisticados, basados en tecnicas Iterativas 
de optimization local, por lo general necesitan de uno o varios parametros de 
regularizacion cuyo valor optimo depende precisamente de la imagen que se 
pretende reconstruir (y que no se conoce), ademas de que la regularizacion 
introducida tiene el efecto de suavizar los contomos de la imagen, haciendola mas 
difusa. Tal es el caso de los metodos que se describen en la patente britanica GB 
2 329 476 del 24 de marzo de 1999. 

Existe pues la necesidad urgente de contar con metodos de reconstruccion 
mejores y con un menor error. Como un ejemplo de ello, recientemente fue 
patentado un metodo de reconstrucci6n de imagenes basado en la aplicacion de 
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10 



redes neuronales artificiales (patente norteamericana US 6,577,700 B1 del 10 de 
junio de 2003). 

En la presente invention, se describen nuevos metodos de reconstruction de 
imagenes basados en recocido simulado y algoritmos geneticos. De acuerdo con 
las leyes de la fisica, en especial de la electrostatica, el sensor usado en TCE (tal 
como se observa en la figura 4) puede ser considerado como un caso particular 
de un sistema de conductores en un medio dielectrico, cuya teoria fue 
desarrollada originalmente por J. C. Maxwell (1873). En nuestro caso, los 
electrodos detectores (B) actuan como los conductores mientras que el tubo 
aislante (A) y el contenido del sensor actuan como el medio dielectrico. Para un 
sensor de n electrodos, la carga electrica q t inducida en los electrodos esta 
relacionadas con el potential electrico v, de los mismos por medio del sigulente 
conjunto de ecuaciones lineales simultaneas: 



15 
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donde los en son los coeficientes de capacitancia propia (o simplemente la 
capacitancia propia) del electrodo i, mientras los demas, eg , con i * j\ son los 
coeficientes de capacitancia mutua (o simplemente las capacitancias mutuas) 
20 entre los electrodos i y jr. Puesta en forma de matrices, la ecuacion (1) queda 
como: 
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Dado que las capacitancias mutuas tienen la propiedad de reciprocidad, esto 
es, cy = Cji, s6lo hay m = in(w-l) capacitancias mutuas independientes, 
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correspondientes a la matriz triangular inferior o superior de C, y correspondientes 
tambien a cada uno de los m pares distintos de electrodos que es posible formar 
en el sensor. 

El valor de las capacitancias mutuas es una funcion no lineal muy 
compticada de la geometria del sistema de conductores y de la distribution 
espacial de la constante dielectrica o permitividad relativa e (que llamaremos 
tambien simplemente 'permitividad') de los materiales dielectricos que los 
separan. En el caso del sensor de TCE, la geometria de los electrodos, la del tubo 
aislante y el valor de su constante dielectrica, son todos fijos. Por lo tanto, 
podemos decir que las capacitancias mutuas son una funcion no lineal s6lo de la 
distribucion espacial de la constante dielectrica en el interior del sensor, eOy). Al 
problema de calcular las capacitancias mutuas correspondientes a una 
distribucion especifica de permitividad dentro del sensor se le llama problema 
directo. 

El empleo de los electrodos de guarda ciHndricos (D) mostrados en la figura 
3, junto con la suposicion de que la distribucion espacial de la permitividad cambia 
poco en el sentido axial, permite representar el sensor por medio de un modelo 
bidimensional (Xie et a/., 1989). A menos que se indique explfcitamente otra cosa, 
en lo sucesivo utilizaremos el modelo bidimensional del sensor. Por lo tanto, las 
cargas electricas q t y las capacitancias c« y Cy de las ecuaciones (1) y (2) seran 
ahora cantidades por unidad de longitud de los electrodos en sentido axial 
Utilizaremos la tilde (o sea, q n c iX y Cy) para denotar las cantidades totales que 

resultan al considerar la longitud real de los electrodos en la practica. Las 
variables anteriores se relacionan entre sf por medio de la longitud de los 
electrodos L, de acuerdo con 



Si dividimos el interior del sensor bidimensional en p regiones de igual area 
(que llamaremos pixeles) donde la permitividad se considera constante, entonces 
la version discreta del problema directo es 



9i = 
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(5) 



donde c es el vector de capacitancias mutuas (por unidad de longitud), fx son 
funciones no lineales que no se conocen explicitamente y z es el vector de 
5 permitividades correspondiente a las p regiones o pixeles dentro de la zona de 
exploration. 

Aplicando la Ley de Gauss, las capacitancias mutuas por unidad de longitud 
axial de los electrodos se pueden calcular como 
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donde e 0 es una constante ffsica llamada permitividad del vacfo, igual a 
8.854 x io 12 faradios por metro, T t es una curva cerrada que rodea al electrodo *, 
dl es un vector (normal) que representa un elemento infinitesimal de la curva I\, 

15 dl es un elemento de longitud de dicha curva, el simbolo 1 • 1 representa el producto 
escalar de dos vectores, y ft es el potencial electrostatico que se produce en el 
sensor al aplicar un voltaje de V voltios en el electrodo j (al que llamaremos 
electrodo fuente o emisor) y 0 voltios en todos los demas (a los que llamaremos 
electrodos receptores). 

20 El potencial <j? esta determinado por la solucion de la siguiente ecuacion 

diferencial parcial 

V-e(x 9 y)V0 J =0 (7) 

25 sujeta a las condiciones de frontera (a) V voltios en el electrodo fuente y (b) 
ft= 0 en los electrodos receptores y en la pantalla externa. En general, la 
ecuacion (7) no tiene solucion analltica y es necesario resolverla numericamente. 

Al problema de estimar cual es la distribucion espacial de permitividad dentro 
del sensor que corresponde a un conjunto determinado de valores de 

30 capacitancia mutua, se le conoce como problema inverso, y es el problema que 
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deben resolver los m§todos de reconstruction de imagenes. Normalmente, la 
estimation de la permitividad se hace en forma discreta, representendola por un 
vector e como el de la ecuacion (5), que se debe calcular a partir de un vector de 
capacitancias mutuas observadas c, obtenido por medio de un aparato de 
medicion apropiado. 

Para resolver el problema inverso, la mayoria de los sistemas de TCE 
emplean el metodo de proyeccion inversa lineal (LBP) (Plaskowski et a/., 1995; 
Yang y Peng, 2003; Xie et a/., 1989, 1992), que se describe a continuation. Como 
primer paso, para cada uno de los m = pares de electrodos que es posible 

formar, se calcula un mapa de sensibilidad definido por 

sm = c ' (fc) " c 'W para i = 1, . . , n (8) 

°i(Ml) C i(emp) 

donde k es el numero de pixel (de 1 a p), Ci(k) es la capacitancia medida con el par 
de electrodos i cuando el area correspond iente al pixel k esta llena de un material 
de alta permitividad y el resto del sensor contiene material de baja permitividad, 
mientras que c^ M //> y c i{ em P ) son las capacitancias del par de electrodos i cuando el 
interior del sensor esta lleno de material de alta y baja permitividad, 
respectivamente. Generalmente estos mapas de sensibilidad se calculan 
resolviendo numericamente la ecuacion (7) y aplicando la ecuacion (6). 

Ya teniendo los mapas de sensibilidad, 6stos se pueden usar para obtener 
una imagen de permitividad a partir de cualquier vector de m = kn(n-X) 
capacitancias mutuas medidas, c, correspondientes a alguna distribution 
desconocida de material dentro del sensor. Para ello se deben primero normalizar 
las lecturas de capacitancia medidas de acuerdo con 

= c '~^> (9) 

C i(/W0 ~~ C i{emp) 

donde A* es la capacitancia mutua normalizada correspondiente al par de 
electrodos i y a es el valor de capacitancia medido en la practica para ese mismo 
par de electrodos. 
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La formula basica del metodo LBP calcula un 'nivel de gris' g(k) para cada 
pixel como 

m 

s(k) = Para k=\,.. ,p (10) 

£*,(*) 

1=1 

En principio, se supone que ese nivel de gris esta relacionado linealmente 
con la permitividad, con g = 1 yg = 0 correspondiendo a la permitividad de los 
materiales de alta y baja permitividad, respectivamente. El metodo LBP se basa 
en realizar una aproximacion lineal de un problema que, como ya se ha 
mencionado, en realidad es esencialmente no lineal (Gamio y Ortiz-Aleman, 
2003). Por lo tanto, este metodo de reconstruction de imagenes provoca errores 
considerables, los cuales son particularmente graves si las diferencias de 
permitividad son grandes. 

Hasta ahora, la principal alternativa al metodo LBP ha sido el empleo de 
metodos iterativos que buscan minimizar una funcion objetivo, usando tecnicas de 
optimization local como el metodo de Newton-Raphson regularizado o variantes 
del mismo (Yang y Peng, 2003). Como un ejemplo de estos metodos esta el que 
se usa en el programa de computo EIT2D (Vauhkonen et aL, 2001), desarrollado 
por investigadores de Finlandia y el Reino Unido. Su metodo de basa en 
minimizar, con respecto a s, el funcional 

Hc meas ~c ca/C || 2 + a 2 I|Ls|| 2 (11) 

donde a es un parametro de regularizacion, L es una matriz que contiene algun 
tipo de informacion a priori sobre la suavidad de c, y c C aic= f(e) es el vector de las 
capacitancias mutuas calculadas para una distribucion de permitividad especifica 
dentro del sensor. En el programa EIT2D, c ca i c se calcula resolviendo el problema 
directo por medio del metodo de elementos finitos. Partiendo de una estimacion 
inicial e G , la minimizacion se realiza por medio del siguiente procedimiento ilerativo 
(basicamente un metodo de tipo Newton con regularizacion de Tikhonov) 

- b* + [J/J* + « 2 L T L] _1 {J/ [c meas - f(s k )]- «VLs t } (12) 
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donde J k se conoce como la matriz jacobiana, Que contiene las derivadas 
parciales de f(s), evaluadas en s k 



J k = J(s,) = 



d£j 



( j=\...p) (13) 



5 Sin embargo, estos metodos de reconstruccion de imagenes tienen el 

problema de que requieren, para su buena operacion, de uno o mas parametros 
de regularizacion cuyo valor correcto tiene una fuerte dependencia, precisamente, 
de las caracteristicas de la distribution de permitividad que se pretende 
reconstruir, lo cual implicaria conocer de antemano la soluci6n del problema a 

10 resolver. Asimismo, generalmente estos metodos producen imagenes 
distorsionadas, debido a que la regularizacion tiene un efecto de suavizamiento 
excesivo sobre la permitividad obtenida. Si la regularizacion es muy fuerte ocurrira 
el suavizamiento mencionado, y si es muy debit el metodo puede volverse 
inestables y/o no converger a la solucion buscada. 

is Estos algoritmos de optimization local, durante su busqueda, exploran solo 

un pequeno sector del dominio de soluciones, restringido a la vecindad que 
circunda la solucion inicial. Si la solucion optima del problema, es decir, el mmimo 
absoluto de la funcion objetivo, se encuentra alejada de la solucion inicial, 
diffcilmente sera alcanzada debido a la presencia de mfnimos relativos 

20 interpuestos en su camino, lugares donde suelen quedar entrampados estos 
metodos. Los metodos mas utilizados correspondientes a esta categoria son la 
inversion lineal por mmimos cuadrados y las tecnicas que emplean el gradiente de 
la funcion objetivo, como el de maxima pendiente (en ingles 'steepest descent') y 
el gradiente conjugado. En general, los metodos de busqueda local explotan la 

25 escasa informacion derivada de la comparacion de una pequena cantidad de 
modelos, evitando asf una busqueda extensiva en el espacio de modelos 
(Sambridge y Drijkoningen, 1992). 

Los metodos de busqueda global, como su nombre lo indica, exploran todo 
el dominio de soluciones a lo largo del proceso de inversion. Hacen un rastreo 

30 exhaustivo en el espacio de modelos. De esta manera, a pesar de existir 
soluciones parciales del problema, es mayor la probabilidad de que la solucion 
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final corresponds al mejor ajuste entre los datos observados y los datos sinteticos. 
Esta clase de metodos, en contraste con las tecnicas locales, no requiere de la 
informacion proporcionada por las derivadas de la funcion objetivo debido a que 
en ellos el problema no requiere ser linealizado. Los algoritmos de optimizacion 
global emplean criterios estocasticos para explorar simultaneamente todo el 
espacio de soluciones en la busqueda del modelo optimo. El mas conocido de los 
metodos globales es el de Monte Carlo, que realiza una busqueda puramente 
aleatoria y no sesgada. Es decir, al generar cada nuevo modelo, no aprovecha la 
informacion obtenida de los modelos previamente evaluados (Gallagher et al., 
1991). La aleatoriedad no encauzada es el rasgo mas caracteristico de este 
metodo, que lo distingue del resto de los metodos globales. Entre las tecnicas de 
optimizacion global, se encuentran tambien los metodos de algoritmos geneticos y 
de recocido simulado. Ambos fueron concebidos como analogias con sistemas de 
optimizacion que ocurren en la naturaleza. Los algoritmos geneticos emulan los 
mecanismos de la evolucion biologica mientras que la tecnica de recocido 
simulado tiene su base en la termodinamica. Ambos metodos son intrinsecamente 
no lineales y, por lo tanto, se prestan naturalmente a su aplicacion a la tomografia 
de capacitancia, un problema no lineal. 

ESPECIF1CACION DE LA INVENCION 

La presente invencion se refiere a un metodo de reconstruction de imagenes 
para la visualization de flujos multifasicos usando tomografia de capacitancia 
electrica, el cual esta basado en tecnicas heurfsticas no lineales de optimizacion 
global, en particular al metodo de recocido simulado y los algoritmos geneticos. 
Asimismo, considera el arreglo circular de electrodos metalicos rectangulares 
contiguos alrededor de la pared externa de un tubo hecho de un material 
electricamente aislante, formando un sensor. A traves este sensor fluye una 
mezcla de fluidos en forma de flujo multifasico, cuya distribucion espacial en el 
interior del sensor se desea conocer. Se obtienen datos a base de efectuar 
mediciones de capacitancia mutua entre todos los electrodos de pared que es 
posible formar. Dichos datos dependen de la distribucion de los fluidos dentro de 
la tuberia o tanque. 
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Por lo tanto, un objeto de la presente invencion es proporcionar un metodo 
para procesar datos a base de efectuar mediciones de capacitancia mutua entre 
todos los electrodos de pared para reconstruir una imagen de la distribution 
espacial de las fases o componentes de la mezcla multifasica, que fluye a traves 
del sensor, utilizando los metodos de recocido simulado y/o de algoritmos 
geneticos. 

Otro objeto de la presente invencion es que se utiliza en medidores 
tomograficos de flujo multifasico para cuantificar los diferentes fluidos (como gas, 
aceite y agua) producidos en un pozo petrolifero. 

Aun otro objeto mas de ia presente invencion es que se usa para la 
visualization y rnonitoreo de otros flujos y procesos multifasicos que ocurren en la 
industria, tales como, optimizacion del diseno y operation de separadores, 
optimizacion del diseno de unidades de desintegracion de fluidizacion catalrdca 
(FCC) y sistemas de lecho fluidizado, y optimizacion del diseno de reactores de 
lecho fijo. 
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BREVE DESCRIPCION DE LOS D1BUJOS DE LA INVENCION 

La figura 1 representa esquematicamente un sistema de tomografia de 
procesos utilizado en la presente invention. 

La figura 2 representa esquematicamente un sistema de tomografia de 
capacitancia electrica empleado en la presente invention. 

La figura 3 muestra un dibujo del sensor empleado en tomografia de 
capacitancia usado en la presente invention. 

La figura 4 muestra un corte transversal del sensor a la altura de la zona de 
medicion utilizado en la presente invention. 

La figura 5 presenta un diagrama esquematico del metodo de recotido 
simulado empleado en la presente invention. 

La figura 6 muestra el diagrama de flujo de la realization del metodo de 
recotido simulado usado en la presente invention. 

La figura 7 presenta un diagrama esquematico del metodo de algoritmos 
geneticos utilizado en la presente invention. 

La figura 8 muestra esquematicamente el proceso de cruza de modelos 
usado en el metodo de algoritmos geneticos empleado en la presente invention. 

La figura 9 muestra esquematicamente el proceso de mutation de modelos 
usado en el metodo de algoritmos geneticos en la presente invention. 

La figura 1 0 muestra el diagrama de flujo de la realization del metodo de 
algoritmos geneticos utilizado en la presente invention. 

La figura 1 1 muestra los resultados obtenidos al reconstruir imagenes de un 
flujo anular simulado de gas y aceite, usando el metodo de recotido simulado de 
la presente invention. 

La figura 12 muestra los resultados obtenidos al reconstruir imagenes de un 
flujo estratificado simulado de gas y aceite, usando el metodo de recotido 
simulado de la presente invention. 

La figura 13 muestra los resultados obtenidos al reconstruir imagenes de un 
flujo simulado de aceite con burbujas de gas, usando el m6todo de recotido 
simulado de la presente invention. 
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DESCRIPCION DETALLADA DE LA INVENCION 

Los nuevos procedimientos de reconstruction de imagenes de la presente 
invention estan basados en metodos heuristicos de optimization global, 
especfficamente el metodo de recocido simulado (llamado en ingles 'simulated 
annealing') y los algoritmos geneticos. 

Una pluralidad de n electrodos metalicos se coloca alrededor de la periferia 
de una region a visualizar y se efectuan mediciones electricas de resistencia o de 
capacitancia entre ellos. O sea, se obtienen m = \n{n-X) mediciones o datos. Las 
mediciones electricas son preferiblemente valores de capacitancia, aunque 
tambien podrian ser valores de resistencia. Como se muestra en la figura 3, los 
electrodos (B) pueden ser de forma rectangular y estar colocados sobre la pared 
externa de un tubo (A) hecho de material electricamente aislante, formando asi un 
sensor, el cual contiene un flujo multifasico o multicomponente (E). El objetivo es 
utilizar dichas mediciones para reconstruir una imagen de la distribution espacial 
del valor de la constante dielectrica (o permitividad) en la region a visualizar. 
Dicha distribution de la permitividad refleja la distribution de las sustancias que 
ocupan el interior del sensor, donde se encuentra la region a visualizar. 
Consideraremos que dicha region esta dividida en p partes iguales. 

Metodo de Recocido Simulado 

El metodo de recocido simulado esta basado en una analogfa con el proceso 
termodinamico de la cristalizacion. Un fluido mineral que se enfria lentamente 
hasta alcanzar un estado de energfa bajo, da lugar a la formation de cristales 
bien definidos. Si, por el contrario, la sustancia abandona su estado de equilibrio 
termico con un enfriamiento subito o partial, entonces el crista! resultante tendra 
muchos defectos, en caso de que la sustancia no forme un vidrio, caracterizado 
por el desorden metaestable de sus moleculas. Este concepto es utilizado en el 
contexto de los m6todos de optimization para reconocer configuraciones o 
modelos potencialmente utiles. 

Los atomos de cada configuration molecular equivalen a los par&metros del 
modelo en el problema inverso (o sea, la permitividad en los distintos pixeles de 
una imagen). La energfa del sistema para dicha configuration se relaciona con la 
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funcion de costo (o de desajuste) asociada al conjunto de parametros que 
intervienen en el modelo. En nuestro caso para la presente invencion, la energia 
del sistema se asocia con la siguiente norma L2: 

m 

TX C Wmeas-c(i)calcY 

E = 4 = -M— ,.,m) (14) 

donde c(z) wea5 son las m capacitancias mutuas medidas y c(i) ca ic son las calculadas 
al resolver el problema directo para una distribucion de perrnitividad determinada 
e. Partiendo de una distribucion inicial de perrnitividad, el metodo genera una 
gama de configuraciones o combinaciones de parametros considerando una 
cierta temperatura T para el proceso. Para este proposito se emplea el criterio de 
Metropolis et al. (1953), que consiste en desplazar un parametro, en cada 
iteracion, una distancia aleatoria y pequena. Este desplazamiento provoca un 
cambio AE en la energia total del sistema. Si AE es menor o igual a cero, el 
desplazamiento del parametro es aceptado y la configuracion resultante se 
considera como la nueva configuracion actualizada. Cuando existe un incremento 
en la energia del sistema (AE es mayor que cero), la probabilidad de aceptacion o 
rechazo para el desplazamiento se determina como: 

P(AB) = e^ r (15) 

Para saber si es o no admitido el cambio de posicion que implica el aumento 
de la energia del sistema, se elige en forma aleatoria un numero entre cero y uno, 
que se compara con el valor de la probabilidad correspondiente a AE. Si es menor 
dicho numero, se admite el desplazamiento y se considera a la nueva 
configuracion como la actualizada; si es mayor, no se admite y se conserva la 
configuracion que se tenia antes del movimiento. Repitiendo sucesivamente este 
procedimiento se esta simulando el movimiento termico de los atomos de un 
sistema (que se encuentra en equilibrio termico), a una temperatura fija T. Si lo 
que se desea es alcanzar el estado base del sistema, es decir, el estado de 
menor energia y mayor ordenamiento, entonces se debe disminuir muy 
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lentamente la temperature para simular un proceso cuasiestatico. Esto significa 
que, durante el enfriamiento, el sistema debe experimentar una sucesion de 
estados infinitesimalmente alejados del estado de equilibrio termico. 

El metodo de recocido simulado tiene tres componentes basicos 
(Vasudevan, 1991): una funcion de energia o funcion de costo, una funcion de 
orden (criterio de Metropolis) y un parametro que controla la temperatura del 
sistema. El proceso consta de tres ciclos anidados. 

En la figura 5 se muestra un diagrama de funcionamiento para el metodo en 
la presente invencion. El ciclo extemo (3) regula la temperatura del sistema. Cada 
vez que se cumple un ciclo, la temperatura del proceso disminuye al ser 
multiplicada por un factor R T que normalmente es cercano a uno (0<i?r<l)- De 
esta manera se lleva a cabo el enfriamiento lento y gradual que se propone. El 
ciclo intermedio (2) se encarga de actualizar los valores, independientes entre si, 
de una serie de constantes Ki asociadas a cada parametro. Dichas constantes 
determinan el mSximo incremento que podra tener cada parametro al momenta de 
ser perturbado en el ciclo mas interno (1 ) del proceso. Los valores que adoptan 
estas constantes dependen de la cantidad de veces que haya sido aceptado el 
modelo actual al termino de cada secuencia de ciclos internos (1) segun el criterio 
de Metropolis. En el ciclo interno (1) se perturban los valores de los parametros 
empleando los factores K x , definidos en el ciclo intermedio (2). La perturbacion se 
realiza multiplicando a cada parametro por el resultado del producto de su 
correspondiente K t por un valor aleatorio entre menos uno y uno. Posteriormente 
se calcula la respuesta sintetica del modelo actual y se evalua el cambio de 
energia en el sistema asociado a la nueva configuration de parametros. Dicha 
variation de energia corresponde al desajuste entre la curva de datos sintetica y 
la observada o medida. Si el desajuste decrece, entonces la nueva configuracion 
sera aceptada como la actual y perturbada de la misma manera. Si, por el 
contrario, la perturbacion aleatoria produce un crecimiento en el desajuste, 
asociado a un incremento en la energia E del sistema, a dicha configuracion se le 
asigna una probabilidad de aceptacion de acuerdo con el criterio de Metropolis. 

Los ciclos (1), (2), y (3) de la figura 5, se repiten, mientras la temperatura del 
proceso disminuye progresivamente. Conforme disminuye la temperatura, las 
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variaciones de los parametros son cada vez mas pequenas. De esta forma, la 
busqueda en el dominio de soluciones tiende a confinarse hacia los modelos 
asociados con el mfnimo absoluto de la funcion de desajuste E. El resultado final 
es un conjunto de valores para los parametros (o sea, la permitividad en los 
distintos pixeles que forman la imagen) cuya respuesta sintetica reproduce los 
datos observados, salvo por un error suficientemente pequeno. 

Como un ejemplo, una realization espetifica, pero no la unica, del metodo 
de recocido simulado se presenta en la figura 6. En el bloque (1), partimos de una 
serie de valores iniciales de la temperatura, las constantes de perturbation, las 
permitividades y la funcion de costo {T 0i K iio) y s m (con i = 1, p) y E 0 ). Al 
mismo tiempo, tambien en el bloque (1) se inicia el contador del ciclo externo, 
denotado por k. Posteriormente, en el bloque (2) se inician los contadores de los 
ciclos interno e intermedio, i y j % respectivamente, y entramos al ciclo intemo. 
Como parte de dicho ciclo, en el bloque (3) se perturba aleatoriamente cada uno 
de los parametros (o permitividades), uno tras otro. Tambien en el bloque (3), 
cada vez que se perturba un parametro, se resuelve el problema directo y se 
calcula la funcion de costo o desajuste, E, aplicando la ecuacion (14). Si hubo un 
decremento de E con respecto a la evaluation anterior, se acepta el valor 
perturbado del parametro como el nuevo valor del mismo, se incrementa el 
contador de parametros i del ciclo interno, y se procede a perturbar el siguiente 
parametro (si lo hay). Si, por el contrario, hubo un incremento de E, entonces se 
aplica el criterio de Metropolis, en el bloque (4), para decidir si se acepta o no el 
valor perturbado del parametro como su nuevo valor actualizado. Si, de acuerdo 
con dicho criterio, se acepta el nuevo valor, entonces se incrementa el contador 
del ciclo interno i y se procede a perturbar el siguiente parametro (si lo hay). Si, de 
acuerdo con el criterio de Metropolis, no se acepta el valor perturbado, entonces 
no se incrementa el contador i y se procede a perturbar nuevamente el mismo 
parametro. Una vez que se han actualizado todos los parametros en el bloque 
(5) se ajusta el valor de las constantes K t (que determinan la forma en la cual se 
perturba a los parametros en el ciclo interno) y se incrementa el contador del ciclo 
intermedio,/. N K determina el numero de veces que se repetira el ciclo interno sin 
que haya una disminucion de la temperatura. O sea, el ciclo intermedio consiste 
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en la repeticion del ciclo interno Nk veces, pero con diferentes valores deK t . Esto 
se hace con el fin de evitar que la temperatura descienda demasiado 
rapidamente, lo cual puede ser contraproducente en algunas aplicaciones del 
recocido simulado. Sin embargo, en la presente invencion en particular no ocurre 
asi, y el valor de N K se toma como uno. Al terminar el ciclo intermedio, en el 
bloque (6) se reduce la temperatura como se indico en parrafos anteriores y se 
incrementa el contador del ciclo externo k. El procedimiento anterior completo se 
repite hasta que k alcanza el Hmite de iteraciones N T , o antes si se alcanza un 
valor suficientemente bajo de la funcion de costo E. 

La presente invencion tambien se refiere a un metodo de reconstruction de 
imagenes basado en algoritmos geneticos, el cual se describe a continuation. 

Metodo de Algoritmos Geneticos 

Los algoritmos geneticos, originalmente propuestos por Holland (1975), 
representan una evolution del metodo de Monte Carlo para problemas de 
optimization fuertemente no lineales. La busqueda del modelo optimo se lleva a 
cabo explorando simultaneamente la totalidad del espacio de soluciones, 
empleando una regla de transition probabilfstica para guiar dicha busqueda. El 
proceso se inicia a partir de un conjunto de modelos selecclonados 
aleatoriamente. 

Los parametros de cada modelo se transforman en codigo binario para 
formar cadenas denominadas cromosomas, sobre las cuales se aplican criterios 
de seleccion natural y genetica. Los procesos de seleccion, cruza y mutacion 
actualizan la poblacion de modelos, dando lugar a una nueva generacion de 
cromosomas, emulando la forma en que los sistemas biol6gicos evolucionan para 
producir prganismos mejor adaptados al entorno. El proceso completo se repite 
hasta que la media de la funcion de ajuste se acerca al maximo ajuste para toda 
la poblacion. 

El diagrama de flujo de la figura 7 resume el proceso utilizado para aplicar un 
esquema de inversion basado en algoritmos geneticos, similar al descrito por 
Rodriguez-Zuniga et al. (1996) y Ortiz-Aleman et al. (2002, 2003) para la inversion 
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de parametros elasticos del subsuelo a partir de datos de inclination del tenreno y 
para la inversion de datos aeromagneticos, respectivamente. Los pasos basicos 
del proceso se describen brevemente a continuation, con referenda a la figura 7. 

Discretization 

En el bloque (1) de la figura 7, los parametros se representan por medio de 
un vector de incognitas s, que representa la distribution de permitividad 
desconocida. La funcion de costo, que determina el ajuste entre los datos 
observados y la respuesta sintetica de un modelo dado, se denota como E(e). La 
codification de los parametros se realiza tomando en cuenta la extension 
necesaria de la busqueda en el espacio de modelos y la resolution deseada 
(Stoffa y Sen, 1991). De esta manera, la extension se define para cada parametro 
estableciendo un par de cotas a t y b t , es decir, a t <%< b t . La resolution se 
controla con el intervalo de discretization d t , definido como 

df = QijJhi (16) 

donde N t es la cantidad de posibles valores para el parametro durante el proceso 
(Sambridge y Drijkoningen, 1992). Los modelos permitidos, s, definidos por el 
conjunto de parametros s t , estan restringidos al dominio de valores 

e i = a i + jd t para j = 0,..., N g (17) 

Poblacion inicial 

Tambien en el bloque (1) de la figura 7, a partir de combinaciones aleatorias 
de los parametros, se construye una poblacion inicial de modelos, cuya dimension 
depende del problema particular a resolver. Cada combination se traduce en un 
conjunto de indices enteros, definidos por la ecuacion (17). Estos valores enteros, 
que establecen el valor particular de cada parametro del modelo, posteriormente 
se codifican como cadenas binarias llamadas cromosomas. Estas cadenas estan 
formadas por grupos consecutivos de bits, llamados genes, que representan el 
valor de los diferentes parametros s t . Los modelos que forman las poblaciones en 
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generaciones posteriores se crean a partir de los tres mecanismos evolutivos 
esenciales: la selection, la cruza y la mutation. 



Problema directo y evaluation de la funcion de costo 

El problema directo, en el bloque (2) de la figura 7, consiste en calcular la 
respuesta teorica o sintetica para cada modelo en el proceso Iterative. 
Posteriormente, en el bloque (3) de la figura 7, dicha respuesta sintetica se 
compara con las observaciones o datos mediante alguna medida de semejanza 
conocida como funcion de costo. El criterio utilizado en este trabajo es la normal 
definida anteriomnente en la ecuacion (14), aunque puede haber otros. La norma 
seleccionada para evaluar el ajuste debe ser sensible a la forma y complejidad de 
las curvas observada y calculada. 



Selection 

En el bloque (4) de la figura 7, a partir de una poblacion de Q individuos y de 
sus respectivas funciones de costo E(s k ) (k=l, ...,Q), se determina para cada 
modelo una probabilidad ('acumulada') de selection P(s k ) que dependera de su 
nivel de ajuste. 

Una formula que puede ser usada para determinar la probabilidad 
acumulada de selection es: 

donde E maxi y E pr0 m son las funciones de costo maxima y promedio de la 
generation, respectivamente, y Q es el numero de individuos de la poblacion. 
Como siguiente paso, se utiliza un procedimiento de ruleta sesgada (Goldberg, 
1989) para seleccionar una nueva poblacion de modelos. Se generan Q numeros 
aleatorios r k entre cero y uno. Si P(s k -i) <r k <P(e k ), entonces e* es seleccionado 
para formar parte de la nueva poblacion. Esto implica que en la nueva poblacion 
puede haber modelos , gemelos\ Adicionalmente, se pueden anadir a la poblacion 
de Q modelos, 'clones 1 de los mejores modelos, los cuales no se someteran a la 
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cruza ni a la mutation, con el objeto de asegurar que dichos modelos deseables 
no se pierdan al pasar poresos procesos que son de naturaleza aleatoria. 

Alternativamente, la probabilidad de selection puede determinarse como 
(Sambridge y Drijkoningen, 1992) 

PCs,) = a-bE(s k ) (19) 
que describe una distribution de probabilidades lineal, y 

P(e k ) = Ae-" M (20) 

que corresponde a una distribuci6n exponential. Los valores que suelen tomar las 
constantes a,b,AyB son los siguientes 

b = Q- X {E max -E prom r , a > bE ntax 

donde E max , E prom y E a son las funciones de costo maxima, promedio, y la 
desviacion estandar de todos los desajustes de la poblacion initial, 
respectivamente. 

Stoffa y Sen (1991) proponen un criterio de selection basado en una 
probabilidad de actualization. El criterio consiste en comparar el desajuste de 
cada modelo perteneciente a la generation actual con el de un modelo de la 
generation anterior, elegido al azar. Si el desajuste del modelo actual es menor, 
entonces este se conserva. De lo contrario, se considera un valor P u que 
establece la probabilidad de sustitucion del modelo anterior por el actual. Este 
procedimiento controla la influencia de los ajustes de generaciones previas sobre 
la poblacion actual. El valor sugerido por Stoffa y Sen (1991) para la probabilidad 
de sustitucion P u es del 90%. 

Cruza 

Los nuevos modelos se engendran en el bloque (5) de la figura 7 a partir de 
una generation progenitora de Q modelos. En forma aleatoria se integran Q/2 
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parejas de individuos. Cada pareja es potencialmente capaz de cruzarse. Para 
determinar cuales de ellas llevaran a cabo la cruza, se le asocia a cada una un 
numero aleatorio entre cero y uno. Si dicho nCimero es menor que la probabilidad 
de cruza P c entonces la pareja correspond iente efectua el apareamiento. De lo 
contrario, la pareja se preserva intacta en la siguiente generation. 

El mecanismo de cruza se basa en elegir aleatoriamente la position de un 
gen para ambas cadenas. Las cadenas se parten en ese punto para intercambiar 
information entre ellas, tal como se observa en la figura 8, donde se muestra en 
(A) la partition de los cromosomas padres en cuatro gametos, y en (B) el 
intercambio y union de dos gametos formando dos cigotos hijos. El proposito de 
cruzar dos cadenas diferentes es explorar nuevas regiones del dominio de 
soluciones, donde pudiera ubicarse el rmnimo absoluto. En el proceso normal de 
cruza, las parejas que se aparean tienen dos hijos, y la poblacion de modelos se 
mantiene automaticamente en Q individuos. Alternativamente, se puede 
considerar que las parejas tienen solo un hijo, en cuyo caso la poblacion se debe 
completar con modelos generados aleatoriamente hasta volver a tener Q 
individuos. 

Mutacion 

La mutacion, al igual que la reproduction sexual (o cruza), propicia la 
diversidad genetica en una poblacion. La mutacion hace posible que la busqueda 
prospere cuando se encuentra confinada en los alrededores de un minimo local. 
La mutacion se realiza en el bloque (6) de la figura 7, mediante el cambio de 
paridad de un bit seleccionado al azar dentro de la cadena binaria (cromosoma). 
El porcentaje de modelos a los cuales se aplica el proceso de mutacion, al igual 
que en la cruza, depende de un parametro denominado probabilidad de mutacion, 
simbolizado por P m . Este mecanismo previene la convergencia prematura del 
metodo, cuando la poblacion es excesivamente homogenea e incapaz de 
continuar el proceso evolutive. En la figura 9 se representa una cadena arbitraria y 
se ejemplifica la mutacion del valor de un bit aleatorio, en (A) se muestra el 
cromosoma original y en (B) el cromosoma mutado. 



25 



WO 2005/019779 



PCT7MX2003/000067 



Una alternativa para definir la probabilidad de mutaci6n P m fue propuesta por 
Yamanaka e Ishida (1996). Consiste en determinar el nivel de homogeneidad de 
los individuos en cada generacion mediante el calculo de un coeficiente de 
variacion promedio 7, para cada parametro, mediante la expresion 



donde p es la cantidad de parametros, s t es el promedio del z-esimo parametro, y 

at es la desviacion estandar. A continuation se define P m como funcion de 7, es 
decir 



donde P- mi es la probabilidad de mutacion initial. 

Con la mutacion concluye la secuencia de operaciones que define a un 
algoritmo genetico. Dicha secuencia se repite hasta satisfacer alguna tolerancia 
preestablecida. 

Como ejemplo, una realization especifica, pero no la unica, del metodo de 
algoritmos geneticos se presenta en la figura 10. En el bloque (1), se parte de una 
poblacion de Q modelos generados aleatoriamente, cada uno representando una 
distribution de permitividad. Al mismo tiempo, tambien en el bloque (1) se inicia el 
contador de iteraciones L A seguir, se inicia el contador de modelos k y comienza 
un ciclo por medio del cual, para cada uno de los Q modelos, en el bloque (2) se 
resuelve el problema directo y se calculan la funcion de desajuste E(e*) y la 
probabilidad acumulada de seleccion P(s k ). Posteriormente, en el bloque (3) se 
realiza la seleccion de modelos con menor E, de acuerdo con el valor de P(s k ) y 
aplicando el procedimiento de ruleta sesgada ya descrito con anterioridad. Luego, 
en el bloque (4) se pasa a la cruza de modelos. En dicho bloque (4), se forman 
aleatoriamente Q/2 parejas, que se cruzan selectivamente de acuerdo con el valor 
de P c , produciendo solo un hijo por pareja. En ese mismo bloque (4), las parejas 
no cruzadas se conservan (ambos integrantes) iguales y la poblacion se completa 




(21) 



P itti para 7>0.1 
-0.1 para 0.02 < 7 < 0.1 
0.2 para y < 0.02 



(22) 
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a Q individuos con nuevos modelos aleatorios. Por ultimo, en el bloque (5) se 
calculan el coeficiente de variaci6n promedio y y la probabilidad de mutacion P m , 
y se realiza la mutacion de modelos de acuerdo con la description dada en 
parrafos anteriores. Todo el procedirniento anterior se repite continuamente, 
incrementando i en cada iteration. El algoritmo se termina cuando se alcanza un 
numero preestablecido de iteraciones N o cuando la funcion de desajuste es lo 
suficientemente pequena. 

Formulation del problema directo. 

El problema directo consiste en el calculo de las capacitancias mutuas c ij9 i* 
j t que resultan de la presencia de una distribution de permitividad 8 en el interior 
del sensor. Tanto el metodo de recocido simulado como el de algoritmos 
geneticos requieren de la solution repetida del problema directo. Por lo anterior, 
es importante contar con un metodo adecuado para resolver dicho problema, que 
logre un equilibrio razonable entre exactitud (o precision) y rapidez. En el contexto 
de esta invention, el problema directo se puede resolver utilizando una rutina 
optimizada desarrollada por los autores basada en el metodo de volumenes 
finitos, que es descrita muy brevemente a continuation. Esta rutina es mas 
eficiente que las reportadas hasta hoy en la literatura (Yang y Peng, 2003) pues 
es comparable en precision con implantaciones basadas en el metodo de 
elementos finitos con mallas de 9,000 elementos triangulares. La velocidad de 
ejecucion es superior a las de los metodos de elementos finitos y diferencias 
finitas. La rutina esta escrita en Fortran 90 y es totalmente transportable (se ha 
probado en computadoras tipo PC y •clusters' de PCs basados en Windows y 
Linux, en estaciones de trabajo tipo SUN y Alpha y en supercomputadoras Cray). 
La rutina puede extenderse al caso tridimensional sin mayores modificaciones y 
es facilmente paralelizable. 

El problema directo se resuelve mediante el metodo de volumenes finitos en 
una configuration cilfndrica. De esta manera se eliminan las soluciones 
indeterminadas en el centro del disco (como sucede en el metodo de diferencias 
finitas) y el refinamiento de la malla se vuelve mas flexible en comparacion con los 
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metodos de elementos finitos. Se resuelven la ecuacion (7), que se repite a 
continuation 

V-e(x,y)V<f> k = 0 

donde e es la permitividad y <j> k es el potential electrostatico que se genera 
cuando el electrodo k es el emisor (o fuente). La ecuacion esta sujeta a las 
condiciones de frontera (a) <j> k =V voltios en el electrodo fuente y (b) $ k = 0 en los 
electrodos receptores y en la pantalla externa. 

Definiendo las coordenadas radial y angular como r y 0, y utilizando el 
metodo de volumenes finitos, la ecuacion discreta es formulada en forma 
conservadora en cada celda % como 

j^<eV^ k )da g = 0 para 1 = 1,..,^ y j=l,..,N e (23) 

donde / y j se refieren a la discretizacion enryft respectivamente, y N r yN e son 
el numero de secciones en las que se subdivide el radio y la circunferencia del 
sensor, respectivamente. 

Aplicando el teorema de Gauss en coordenadas polares, las ecuaciones 
discretas pueden escribirse como 

J r eV^-dTy = 0 (24) 

donde es la frontera de la celda de volumen finito £\ . La frontera Ty se define 
por T w y T E a lo largo de las coordenadas radiales, y y r 5 a lo largo de las 
coordenadas angulares. La ecuacion (24) puede expresarse como la suma de los 
flujos a traves de las caras T N , T s , T E y T w . 
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De (25), el termino correspondiente a los flujos en el radio cero desaparece y 
el problema es equivalente a resolver las ecuaciones en la cercania del centro 
sobre triangulos que tienen un vertice en el centro. Entonces, el sistema discreto 
de ecuaciones para el problema directo resulta bien planteado. El sistema 
completo es similar a un sistema de ecuaciones Laplaciano y debe resolverse un 
sistema diagonal bandeado que incluye las condiciones de frontera periodicas 
impuestas por la geometria del problema. La matriz correspondiente es positiva 
definida y no simetrica, caracteristicas que se pueden aprovechar para 
seleccionar especfficamente los metodos optimos de solucion. 

Finalmente, las capacitancias mutuas se calculan integrando los gradientes 
del potencial a lo largo de una curva cerrada que rodea los electrodos, de acuerdo 
con la ecuacion (6), que se repite a continuation 

v k V I V I dn 

donde e G es la permitividad del vacio (8.854 x io~ 12 faradios por metro), T t es una 
curva cerrada que rodea al electrodo i, dl es un vector (normal) que representa un 
elemento infinitesimal de la curva F h dl es un elemento de longitud de dicha curva 
y $ k es el potencial electrostatico que se produce en el sensor al aplicar un voltaje 
de V voltios en el electrodo k (fuente) y cero voltios en todos los demas 
(receptores). La integration se realiza utilizando una regla trapezoidal y los 
gradientes de potencial se calculan al cuarto orden. 

Durante el procedimiento de reconstruction de una imagen de permitividad 
por el metodo de recocido simulado, es necesario resolver el problema directo y 
encontrar el potencial electrico repetidamente para sucesivas distribuciones de 
permitividad relativamente similares, mientras el metodo converge a la solucion 
final. En virtud de que el potencial correspondiente a dichas distribuciones 
sucesivas cambia relativamente poco, es posible acelerar todo el proceso si se 
emplea un metodo iterativo para resolver el problema directo, tomando como 
estimacion inicial del potencial a la solucion correspondiente a la configuracion de 
permitividad anterior. Como la estimacion inicial del potencial estara muy cerca de 
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la solucion, dicho metodo iterative convergera en pocas iteraciones, alcanzando 
rapidamente una exactitud aceptable. 

EJEWIPLO 

Con objeto de evaluar el desempeno de los metodos de reconstruction de 
imagenes descritos anteriormente, calculamos un conjunto de datos TCE 
sinteticos con la subrutina del problema directo. Para emular las fuentes 
principales de incertidumbre en los datos producidos por uri instrumento de 
medicion operando en condiciones de trabajo (esto es, los errores aleatorios 
inherentes al proceso de medicion de cualquier variable fisica y los errores 
debidos a la precision propia del instrumento de medicion), los calculos de las 
capacitancias sinteticas para el modelo ideal fueron evaluadas con una precision 
numerica del orden de 10" 11 en el metodo iterativo utilizado en el problema directo 
para el calculo de los potenciales. Con objeto de simular la imprecision 
(sistematica) asociada al sensor TCE, durante el proceso de inversion (estimation 
de la distribution de permitividad electrica en el interior del tubo) se utilizo una 
precision sensiblemente menor (KT 5 ) en el calculo de los potenciales en el 
problema directo. 

Las interpretaciones sin restricciones de datos de campos potenciales 
pueden ser de muy poco interes practico debido a la gran ambiguedad presente 
entre las observaciones y las soluciones estimadas. La no unicidad en los 
problemas con campos potenciales surge principalmente de dos fuentes: la 
primera es la ambiguedad inherente provocada por la fisica del problema que 
permite la existencia de una infinidad de soluciones que reproducen la anomalia 
de campo potential; la segunda resulta de la utilization de un numero finite de 
datos que estan contaminados por errores y que pueden contener informaci6n 
insuficiente para construir una solucion unica del problema. Las estrategias que 
permiten superar esta no unicidad consisten en la incorporation de suficiente 
informacion a priori para restringir las soluciones resultantes a una region del 
espacio de parametros que es considerada fisicamente razonable (Pilkington, 
1997). En el caso particular de la tomografia de capacitancia electrica se cuenta 
con informacion acerca de los valores tfpicos de la permitividad electrica para el 
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gas, aceite y agua. El conocimiento de estos valores, asf como el contraste 
significativo entre las propiedades del gas y el aceite, permiten aplicar las tecnicas 
de inversion que aqui se discuten a datos de TCE para flujos bifasicos 
contaminados por errores relativos menores al 2% (siendo este el nivel maximo 
de error que producen los sistemas de adquisicion de datos empleados 
normalmente en TCE). Por otra parte, si se cuenta con una estimation estadistica 
precisa de las incertidumbres en los datos, entonces es posible construir un 
modelo que considere las incertidumbres en los datos y en las estimaciones de 
los parametros, utilizando el esquema que proponen diversos autores para la 
inversion de datos de campos potenciales (por ejemplo, Sen y Stoffa, 1995). 

Usando la subrutina del problema directo a base de volumenes finitos, se 
calcularon las capacitancias sinteticas considerando tres distribuciones 
especificas de permitividad, correspondientes a tres patrones de flujo tipicos 
(anular, estratificado y de burbujas). Simulamos un sensor TCE de 12 electrodosy 
los valores de capacitancia para todas las combinaciones posibles de electrodos 
fueron calculados. Esos fueron los datos simulados. Consideramos una 
distribution de dos componentes con un material de permitividad baja de 1 (gas) y 
otro de permitividad alta de 2.5 (aceite). Nuestras pruebas numericas se han 
hecho con datos libres de ruido y con datos contaminados con errores aleatorios 
de hasta 1%. El algoritmo esta codificado en Fortran 90 y corre en una 
computadora Pentium 4 de 1.7 gigahertz con 512 megabytes de memoria RAM. 
Probamos con una malla de 120*60 elementos para reducir los tiempos de la 
inversion (~ 30 minutos para 60,000 iteraciones), pero los resultados son validos 
para dimensiones mas grandes. La validez de los resultados de estas 
simulaciones fue corroborada con resultados experimentales de laboratorio 
obtenidos empleando modelos fisicos y un sistema real de tomografia de 
capacitancia. 

Despues de una parametrizacion adecuada, tanto el metodo de algoritmos 
geneticos como el de recocido simulado produjeron resultados muy similares, 
para los tres patrones de flujo considerados. La calidad de las imagenes de 
permitividad reconstruidas depende principalmente del numero de iteraciones del 
metodo, tal y como sucede para muchas otras aplicaciones (por ejemplo, Ortiz- 
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Alem£n et a/., 1999, 2001, 2002, 2003; Cruz-Atienza, 1999). En la figura 11 se 
presenta la imagen ideal o Verdadera' de un flujo anular simplrficado (A), y las 
reconstrucciones obtenidas despues de 30,000 (B) y 60,000 (C) calculos del 
problema directo (o iteraciones del metodo). En las figuras 12 y 13 se presentan 
figuras similares para un flujo estratificado y un flujo de burbujas, 
respectivamente, mostrando en (A) la imagen ideal, en (B) la imagen reconstruida 
con 30,000 iteraciones y en (C) la imagen reconstruida con 60,000 iteraciones. En 
todos los casos mostrados en las figuras se uso el metodo de recocido simulado 
(con los algoritmos geneticos no hubo diferencias significativas en los resuttados) 
y se partio de una distribution de permitividad initial homogenea. Estos 
resultados muestran claramente que las imagenes reconstruidas son muy 
parecidas a las imagenes ideales o Verdaderas'. La exactitud de estas 
reconstrucciones es considerablemente mejor que las reportadas hasta ahora en 
la literatura (Yang y Peng, 2003). Si bien los metodos de esta invention no 
requieren, para converger, que la estimation initial de la distribution de 
permitividad este proxima a la solution, es posible acelerar un poco el proceso si 
se emplea como estimation inicial la imagen resultante de un metodo directo 
simple como el LBP. 
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NOVEDAD DE LA INVENCION 

Habiendo descrito la presente invention, esta se considera como novedad y 
por lo tanto, se reclama como propiedad el contenido en las siguientes clausulas: 

1 Un metodo de reconstruction de imagenes que comprende: (a) la obtencion 
de datos medidos correspondientes a la imagen; (b) el procesamiento de 
dichos datos medidos por medio de metodos heuristicos no lineales de 
optimization global, espetificamente el metodo de recocido simulado o los 
algoritmos geneticos, para obtener una imagen; y (c) el desplegado de dicha 
imagen en un dispositivo de visualization. 

2. - Un metodo de conformidad con la clausula 1 , caracterizado porque los datos 

medidos son parametros electricos en el perimetro de una region tal como el 
interior de una tuberia, el de un pozo o el de un tanque, para obtener una 
imagen de la distribution espacial de la permitividad electrica (o conslante 
dielectrica) o de la conductividad electrica en la que se refleja la distribution 
espacial de los materiales o sustancias, tales como gases y/o Hquidos. 

3. - Un metodo de conformidad con la clausula 2, caracterizado porque los 

parametros electricos son valores de capacitancia electrica medidos entre 
los distintos electrodos un sensor formado por una pluralidad de electrodos 
colocados a lo largo del perimetro de la region (tuberia, pozo, tanque). 

4. - Un metodo de conformidad con las clausulas 2 y 3, caracterizado porque la 

imagen obtenida es una distribution espacial de la permitividad electrica (o 
constante dielectrica) en la region de interes (tuberia, pozo, tanque), la cual 
refleja la distribution espacial de los materiales o sustancias, tales como 
gases y/o Hquidos, que ocupan dicha region (tuberia, pozo, tanque). 

5. - Un metodo de conformidad con las clausulas 2 a 4, caracterizado porque la 

imagen esta formada por una cantidad finita de subregiones o pixeles en la 
visualization, cuyo numero depende de la resolution deseada. 

6. - Un metodo de conformidad con la clausula 3, caracterizado porque el sensor 

esta formado por un tubo fabricado de un material electricamente aislante, 
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sobre cuya pared externa se coloca un anillo de electrodos metalicos 
rectangulares. 

7. - Un metodo de conformidad con las clausulas 3 y 6, caracterizado porque el 

sensor contiene un flujo multifasico o multicomponente, y obtiene imagenes 
que muestran la distribution de las fases o componentes, tales como gases 
y/o liquidos de acuerdo a lo dicho en antecedentes. 

8. - Un metodo de conformidad con la clausula 1 , caracterizado porque emplea 

el metodo de recocido simulado como el metodo de optimization. 

9. - Un metodo de conformidad con la clausula 8, caracterizado porque se 

minimiza iterativamente una funtion de costo asociada con la energia del 
sistema, con respecto a e (el vector de los valores de la permitividad en cada 
pixel de la imagen), y siendo de la forma: 

±[cr-cr lc (* k )] 2 

E k = Lm = — " (i=l,..,m) 

i=i 

donde c™ eas son las m capacitancias mutuas medidas y cf c (z k ) son las 
calculadas al resolver el problema directo para una distribution de 
permitividad determinada e*. 

10. - Un metodo de conformidad con la clausula 9, caracterizado porque para 

efectuar la minimization se emplea el criterio de Metropolis. 

11. - Un metodo de conformidad con la clausula 10, caracterizado porque en el 

proceso iterativo de minimization se emplea como estimation initial de la 
distribution de permitividad, 8, tanto una distribution homogenea como la 
que resulta de aplicar a los datos medidos el metodo de proyeccion inversa 
lineal (LBP). 

12. - Un metodo de conformidad con la clausula 9, caracterizado porque el calculo 

de las capacitancias calculadas cf c (e k ) , conocido como el problema directo, 

se realiza por medio del metodo de volumenes finitos. 

13. - Un metodo de conformidad con la clausula 12, caracterizado porque para la 

solution del sistema de ecuaciones que resulta al resolver el problema 
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directo (calculo de cf^s*)), se utilizan metodos iterativos que convergen 

rapidamente al emplear como estimacion inicial del potencial electrostatico, 
el potencial que resulto de la solucion del problema directo en la iteration 
anterior del problema inverso. 

14. - Un m6todo, de conformidad con )a clausula 1, caracterizado porque se 

emplea el metodo de algoritmos geneticos como el metodo de optimization. 

15. - Un metodo de conformidad con la clausula 14, caracterizado porque, a partir 

de una poblacion inicial de Q modelos de permitividad 6^ 0) ...,(?), 
aplicando mecanismos evolutivos como seleccion, cruza y mutation, se 
obtienen iterativamente nuevas poblaciones. 

16. - Un metodo de conformidad con la clausula 15, caracterizado porque los 

individuos e k poseen una funcion de costo o de desajuste pequena, donde la 
funcion de desajuste esta dada por: 

per -cr lc M] 2 

E k = L 2 W = — (1 = 1,.., m) 

/=i 

donde c" ieas son las m capacitancias mutuas medidas y ^(e^) son las 

calculadas al resolver el problema directo para el modelo (o distribucion de 
permitividad) e*. 

17. - Un metodo de acuerdo con las clausulas 15 y 16, caracterizado porque la 

probabilidad acumulada de seleccion para un modelo determinado e* esta 
dada por 



Q ( E max ~Eprom) 



donde E max , y E prom son las funciones de costo maxima y promedio de la 
generacion, respectivamente, y Q es el numero de individuos de la 
poblacion. 
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18. - Un metodo de conformidad con la clausula 17, caracterizado porque se 

puede emplear el procedirniento de la ruleta sesgada para decidir cuales 
modelos son seleccionados en cada iteration. 

19. - Un metodo de conformidad con la clausula 18, caracterizado porque la cruza 

y la mutacion de modelos se realizan aleatoriamente de acuerdo las 
probabilidades de cruza y de mutacion, P c y P m . 

20. - Un metodo de conformidad con la clausula 19, caracterizado porque la 

probabilidad de mutacion P m esta determinada mediante el coeficiente de 
variation promedio y, dado por 



donde p es la cantidad de parametros, s i es el promedio del /-esimo 

parametro, y oi es la desvtacion estandar. 
21.- Un metodo de conformidad con las clausulas 19 y 20, caracterizado porque 
P m se define como funcion de es decir 



donde P ini es la probabilidad de mutacion inicial. 
22.- Un metodo de conformidad con las clausulas 15 y 16, caracterizado porque 
el calculo de ^(e*), conocido como el problema directo, se realiza por 
medio del metodo de volumenes finitos. 




P bli para />0.1 

0.1 para 0.02 <;r< 0.1 

0.2 para y < 0.02 
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